##########################################################################################

library(data.table)
library(dplyr)
library(ggplot2)
library(optparse)
library(clusterProfiler)
library(org.Hs.eg.db)
library(tidyverse)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--diff_exp_file"),type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)


diff_exp_file <- opt$diff_exp_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)
setwd(out_path)

##########################################################################################
# 读取文件
dat_diff <- data.frame(fread(diff_exp_file))

##########################################################################################
## 提取差异高表达基因
result_diff <- dat_diff

computGO <- function(gene_symbol = gene_symbol , type = type){
    gene_df <- bitr(gene_symbol,fromType = 'SYMBOL',toType = 'ENTREZID', 
                    OrgDb = 'org.Hs.eg.db')
    ego <- data.frame(enrichGO(gene_df$ENTREZID,OrgDb = org.Hs.eg.db,ont = 'ALL'))
    ego$type <- type
    return(ego)
}

## 通路富集
enrich_df <- c()
for( clust in unique(result_diff$CLUSTER) ){
    print(clust)
    gene_symbol <- unique(subset( result_diff , CLUSTER == clust )$Hugo_Symbol)
    tmp_go <- computGO(gene_symbol = gene_symbol , type = clust)
    enrich_df <- rbind( enrich_df , tmp_go )
}

write.table(enrich_df, 'enrich.go_filt.txt', sep = '\t', row.names = FALSE, quote = FALSE)

##########################################################################################
## 每种类型的通路取前5最显著的
result <- c()
top_n <- 5
for(t_type in unique(enrich_df$type)){
  for( t_path in unique(enrich_df$ONTOLOGY)){
    tmp <- subset( enrich_df , type == t_type & ONTOLOGY == t_path )
    tmp <- tmp[order(tmp$qvalue , decreasing=F),][1:top_n,]
    result <- rbind(result , tmp)
  }
}

## 只看BP通路
result <- subset(result , ONTOLOGY=="BP")

gene_ratio <- as.numeric(sapply(strsplit(result$GeneRatio , "/") ,"[" , 1))/as.numeric(sapply(strsplit(result$GeneRatio , "/") ,"[" , 2))
bg_ratio <- as.numeric(sapply(strsplit(result$BgRatio , "/") ,"[" , 1))/as.numeric(sapply(strsplit(result$BgRatio , "/") ,"[" , 2))
result$OR <- gene_ratio/bg_ratio
result$y_p <- -log10(result$p.adjust)
result <- result[order(result$y_p , decreasing = T),]
result$Description <- factor( result$Description , levels = unique(result$Description) )

plot <- ggplot(result , mapping = aes(x=Description, y=y_p)) +
    facet_wrap(vars(type) , nrow = 1 , scales = "free_x") +
    geom_bar(stat = 'identity', position = 'dodge') +
    theme_bw()+
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.text = element_text(size = 20),  # 设置分面标签文本大小
        axis.title.x = element_blank(),
        axis.text.x = element_text(color = "black", size = 16, angle = 45, hjust = 1),
        axis.text.y = element_text(color = "black", size = 14),
        axis.title.y = element_text(color = "black", size = 16),
        legend.position = "top",
        legend.title = element_text(color = "black", size = 13),
        legend.text = element_text(color = "black", size = 13))+
          labs(x='GOBP', y='-log10(qvalue)')

out_name <- paste0( "enrichGO.pdf" )
ggsave(filename = out_name , plot = plot, width = 15, height = 10)

